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UNCERTAINTY  PREDICTION  IN 
PASSIVE  TARGET  MOTION  ANALYSIS 

STATEMENT  OF  GOVERNMENT  INTEREST 
[0001]  The  invention  described  herein  may  be  manufactured  and 
used  by  or  for  the  Government  of  the  United  States  of  America 
for  governmental  purposes  without  the  payment  of  any  royalties 
thereon  or  therefor. 

CROSS  REFERENCE  TO  OTHER  PATENT  APPLICATIONS 
[0002]  None . 

BACKGROUND  OF  THE  INVENTION 

(1)  Field  of  the  Invention 

[0003]  The  present  invention  is  directed  to  a  target  motion 
analysis  system  and  method  and  more  particularly  to  a  target 
motion  analysis  system  that  indicates  uncertainty. 

(2)  Description  of  the  Prior  Art 

[0004]  In  the  bearing  only  target  motion  analysis  (TMA) 
problem,  one  must  estimate  the  position  and  velocity  of  a 
constant  velocity  target  from  bearing  observations.  FIG.  1 
shows  use  of  a  target  motion  analysis  system.  A  target  10  is 
positioned  at  an  unknown  location  and  following  an  unknown 
course  relative  to  an  observer  12.  Observer  12  has  a  sensor 
array  such  as  a  passive  sonar  or  radar  array  14  that  can  receive 
a  signal  radiating  from  target  10.  Observer  12  also  has 
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environmental  sensors  16  and  a  database  18.  Array  14,  sensors 
16,  and  database  18  are  all  connected  to  a  processor  20.  User 
input  22  can  also  be  provided  to  processor  20.  Processor  20  is 
capable  of  beamforming  data  from  the  sensor  array  14  in  order  to 
provide  bearings  from  the  observer  12  to  the  target  10. 

Processor  20  can  also  utilize  information  from  the  array  14, 
database  18,  user  input  22  and/or  environmental  sensors  16  to 
provide  an  estimate  of  the  bearing  and  range  to  target  10. 

[0005]  The  target  motion  analysis  system  24  receives 
information  from  processor  20  to  provide  a  TMA  solution.  This 
solution  is  comprised  of  parameters  such  as  the  range,  bearing, 
course  and  speed  of  the  target.  The  target  motion  analysis 
system  24  can  be  a  routine  implemented  on  processor  20.  The 
solution  is  provided  for  display  on  a  suitable  display  device  or 
for  use  by  other  systems. 

[0006]  As  shown  in  FIG.  2,  target  range  R  and  bearing  B 
define  the  position  of  the  target  10  relative  to  that  of  the 
observer  12,  and  target  course  C  and  speed  S  define  the  velocity 
of  the  target  10  as  it  transits  through  the  space.  Target  state 
is  defined  as  x(t)  =  [x(t),y(t),x, y]T  where  x(t)  and y(t)  are  Cartesian 
coordinates  of  the  target  with  respect  to  time,  and  x  and  y  are 
components  of  the  target  velocity. 

[0007]  As  is  well  known,  a  fundamental  property  of  bearings- 
only  target  motion  analysis  (TMA)  is  that  bearing  B  to  the 
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target  10  results  directly  from  sensor  system  output,  but  the 
range  R  to  the  target  10  must  be  inferred  from  additional  data. 
This  additional  data  can  be  obtained  by  maneuvering  the  observer 
12  having  the  sensor  system.  A  single-leg  of  observer  motion 
only  provides  a  partial  target  motion  solution.  (A  leg  is 
defined  as  a  time  interval  of  constant  platform  velocity.) 

[0008]  In  addition  to  a  point  estimate  of  the  current  target 
state,  a  representation  of  the  resulting  uncertainty  in  the 
solution  is  required.  The  positional  uncertainty  is  usually 
defined  by  an  area  of  uncertainty  (AOU)  about  the  target  10 
shown  as  dashed  line  28.  Uncertainty  characterizations  regarding 
other  parameters  may  be  required  as  well.  For  example,  the 
uncertainty  related  to  target  velocity,  comprised  of  the 
target's  course  and  speed,  can  also  be  indicated. 

[0009]  Observations  of  contact  state  are  imperfect.  This  can 
be  modeled  by  representing  the  observation  as  a  mapping  from  the 
true  contact  state  corrupted  by  additive  noise. 

zm  —  ”1”  Vm  (1) 

where  x(tm)  denotes  the  target  state  at  time  tm  from  which  the 

estimated  bearing,  zm,  is  measured  as  a  function  of  the  true 

bearing,  Sm(x(tm)).  The  observation  error,  is  assumed  to  be 

independent  and  identically  distributed  over  the  ensemble  of 

2 

observations ,  >7~N(0,cr),  with  a  denoting  the  standard  deviation  of 
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the  error.  An  estimate  of  the  target  state,  x,  can  be  developed 
by  processing  a  collection  of  observations,  Z  —  \z1,z2,  using 

any  of  a  number  of  commonly  accepted  estimation  methods.  For  an 
estimated  target  track  to  be  considered  plausible,  it  must  fit 
the  data  reasonably  well.  A  standard  measure  of  the  goodness- 
of-fit  is  the  sum  squared  error. 


M 


m  =  ^  Ors 


xm) 


m= 1 


(2) 


[0010]  The  mathematical  relationships  denoting  dynamic 
vehicular  motion  and  the  mapping  between  target  states  and  the 
measurements  used  to  estimate  them  are  often  non-linear.  This  is 
true  for  the  bearing  observation: 


Bm  =  tan  1 


fVx  (tjn)  \ 
yyO-m )  ) 


(3) 


where 


T xO-m )  x(^m)  x0  0-m)  (4) 

TyO-m)  —  y(tjn)  YoC^m)  (3) 

with  (x,  y)  and  (x0,  y0)  denoting  the  coordinates  of  the  target  and 
the  observer  at  time  tm.  Hence,  direct  analytical  expressions  for 
estimated  values  cannot  always  be  obtained  and  iterative 
estimation  processes  must  be  developed  in  order  to  converge  to 
an  acceptable  solution.  One  such  method  applies  a  Gauss-Newton 
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approach  to  the  iteration  process.  For  these  types  of 
algorithms,  the  state  vector  estimated  is  adjusted  with  each 
iteration  by  an  amount  proportional  to  the  change  in  cost 
function,  as  in  equation  2,  that  is  dictated  by  the  gradient  at 
the  current  estimate.  This  iteration  takes  the  form 


xi+i  =xi+s-gi  (6) 

where  i  denotes  the  iteration  number  and  s  denotes  a  scale  factor 
referred  to  as  the  step-size.  When  the  measurement  error  over  a 
batch  of  data  can  be  assumed  to  be  independent  and  identically 
distributed,  a  Gauss-Newton  approximation  to  the  cost  function  g 
will  yield  the  vector  equations: 

gt  —  [ata]~1at(z  —  z)  (7) 

where  A  denotes  an  (M  x  N)  gradient  matrix  of  partial 
derivatives  defining  the  change  in  measurement  with  the  change 
in  state. 
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native  coordinates  onto  a  coordinate  system  with  a  single  non¬ 
zero  component.  When  applied  to  a  matrix,  it  is  only  the  first 
column  that  is  reflected  to  a  single  entry  basis.  Hence,  the 
transformation  is  applied  sequentially  until  each  component  of 
the  state  vector  has  gone  through  a  reflection  process.  The 
combined  transformation  yields  a  transformation  matrix,  T,  with 
the  property  that  TrT  =  I.  the  N-dimensional  identity  matrix. 

This  transformation  matrix  T  is  not  actually  calculated,  only  the 
result  of  the  transformation  is  calculated.  When  applied  to 
gradient  matrix  A,  the  result  is  a  partition  into  a  non-zero 
upper  triangular  square  matrix  Rand  a  remaining  component 
comprised  of  zero  values. 


Assuming  A  is  non-singular  (i.e.,  that  all  the  components  of  x 
are  observable  from  the  data) ,  then  application  of  the 
Householder  transformation  into  equation  7  yields  the  back 
substitution  equation: 

gt  =  [AtTtTA]~1AtTtT(Z  ~Z)  =  R1R  TRTze  =  R~xze  (10) 

The  error  vector  ze  denotes  the  first  n  components  of  the  result 
obtained  by  applying  transformation  T  to  the  residual  vector 
T  (Z  —  Z)  =  [zezr  ]  .  The  upper  triangular  nature  of  the  matrix  R 
means  that  the  components  g,-  can  be  solved  for  sequentially  via  a 
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back  substitution.  Iteration  in  Xj  continues  until  the 
improvement  in  J  (x)  becomes  insignificant  in  the  applied 
termination  criterion. 

[0012]  In  many  state  estimation  problems,  the  representation 
and  visualization  of  uncertainty  in  state  estimates  is  a 
significant  part  of  the  estimation  process.  In  the  case  of 
Kalman  filter  and  its  many  extensions,  the  uncertainty  of  states 
is  assumed  to  be  Gaussian.  While  this  is  useful  computationally, 
it  is  very  often  inadeguate  for  characterizing  the  accuracy  of 
estimates  when  data  quality  is  low  and  the  estimation  problem  is 
nonlinear.  In  recent  years,  alternatives  have  been  developed. 
Unscented  Kalman  filter  have  improved  uncertainty 
characterization.  Particle  filters  have  addressed  this  issue  by 
calculating  the  uncertainty  numerically  through  the  ensemble  of 
states  represented  by  the  particles.  These  particles  can  then  be 
projected  to  lower  dimensional  spaces  to  illustrate  the 
uncertainty  to  an  end  user.  However,  particle  representations  of 
uncertainty  can  be  very  expensive  computationally. 

[0013]  Another  alternative  are  purely  grid  based  methods 
which  consider  a  discrete  subset  of  the  state  space  as  potential 
solution  and  then  produce  error  surfaces  at  the  prescribed 
resolution.  If  well-conceived,  grid  based  approaches  can  allow 
an  efficient  view  of  state  uncertainty,  even  with  some  of  the 
parameters  being  unobservable.  In  the  case  of  bearings-only 
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target  motion  analysis  the  Parameter  Evaluation  Plot  (PEP)  is 
one  example  of  such  a  grid-based  approach.  U.S.  Patent  No. 
7,020,046  discloses  one  version  of  this  method  and  is 
incorporated  herein.  The  PEP  algorithm  utilizes  the  endpoint 
coordinate  system  defining  the  position  of  the  target  relative 
to  the  observer  at  each  of  two  times.  The  base  algorithm  has 
been  a  standard  method  for  target  motion  analysis  (TMA)  for  many 
years  and  various  approaches  have  been  applied  to  improve  it. 

This  is  shown  in  FIG.  3  where  Bi  and  Ri  denote  the  bearing  and 
range  from  observer  12A  to  the  target  10A  at  time  ti,  and  B2  and 
R2  denote  the  bearing  and  range  from  observer  12B  to  target  10B 
at  time  t2  as  indicated  in  FIG.  3.  Between  times  ti  and  t2, 
observer  12A  has  traveled  over  the  path  shown  at  16.  Assuming 
constant  direction  and  velocity  for  target  12A  gives  the 
estimated  path  18  over  the  same  time  period.  For  bearings  only 
estimation  in  endpoint  coordinates,  x  —  [B1,R1,B2,R 2]  must  be 
determined.  The  position  of  the  target  10  relative  to  the 
observer  12  at  arbitrary  time  tm  is  given  by: 

rx(tm)  =  AT1R1  sinCBj)  +  A T2R2  sin(B2)  +  x0(tm)  (11) 

ry(tm )  =  A T±R±  cosCBj)  +  A T2R2  cos(B2)  +  y 0(tm)  (12) 

Here  A Tt  =  (t2  -  tm)/AT  and  A T2  =  (tm  -  tJ/AT  denote  the  time 
displacement  between  the  times  that  the  observations  are  made 
and  the  respective  endpoint  timelines,  normalized  by  the  time 
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interval  between  the  timelines,  A T  =  (t2  —  1 1)  .  The  observer 
displacements : 

Ax0(tm)  =  Ar2x0(t1)  +  ATiXoCtz)  -  x0(tm)  (13) 

Ay o(tm)  =  Ar2y0( +  AT^yoCtz)  -  y0(tm)  (14) 

denote  deviations  of  the  observer  trajectory  from  a  straight 
line  between  the  observer  positions  at  the  respective  time 
lines.  Hence  the  bearing  relation  expressed  in  endpoint 
coordinates  becomes 


_  tan-i  f  A7i^i  sin(Bi)  +  Ar2^2  sin(B2)  +  x0(tm)  \  (15) 

m  311  W.R,  cos(Bi)  +  A T2R2  cos(B2)  +  y0(tm)/' 

[0014]  When  the  endpoint  bearings  Bi  and  B2  are  assumed  to  be 

known  and  fixed  values  can  be  assigned  to  them,  the  PEP 

estimation  problem  reduces  to  two-dimensional  search  for  the 

best  Ri  and  R2  tracking  parameters.  In  this  case,  the  goodness- 

of-fit  can  be  evaluated  directly  over  the  endpoint  range 

parameter  space.  Goodness-of-f it  is  measured  according  to  the 

RMS  performance  criterion  given  by: 


RMS(x )  = 


(16) 


[0015]  This  approach  is  intuitive  in  that  it  considers  fixed 
values  for  directly  observed  states  and  calculates  the  cost 
function  over  a  gridded  sampling  of  the  unobserved  states.  The 
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resultant  cost  function,  shown  in  FIG.  4,  illustrates  a  two- 
dimensional  slice  of  the  higher  dimensional  cost  function.  The 
problem  with  this  approach  is  that  choosing  the  fixed  values  can 
be  complicated  with  no  single  set  of  values  being  representative 
over  the  whole  state  space.  That  is,  the  selection  of  two 
observed  values  is  sensitive  to  the  accuracy  of  the 
observations.  Generalized  statements  can  be  made  about  FIG.  4. 
Due  to  practical  constraints,  extreme  range  variations  between 
ranges  Ri  and  R2  are  unlikely  resulting  in  a  high  cost  function 
near  the  axis.  A  low  cost  function  is  given  by  intermediate  low 
ranges  with  little  variation  since  these  smaller  range  changes 
are  more  plausible. 


SUMMARY  OF  THE  INVENTION 

[0016]  It  is  a  first  object  of  the  present  invention  to 
provide  a  bearings-only  target  motion  analysis  method. 

[0017]  Another  object  is  to  provide  such  a  method  that 
accurately  expresses  the  uncertainty  caused  by  weakly  observable 
parameters . 

[0018]  Yet  another  object  is  to  provide  such  a  method  that 
uses  likelihoods  to  compute  target  motion  and  uncertainty  in  an 
efficient  manner. 

[0019]  Accordingly,  there  is  provided  a  method  for  target 
motion  analysis  for  calculating  the  track  between  an  observer 
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and  a  target  includes  initializing  strongly  observable  state 
components  into  strong  state  vectors  and  weakly  observable  state 
components  to  form  a  weak  state  analysis  grid.  Initial 
probability  density  functions  are  calculated  over  the  weak  state 
grid  using  strong  state  components.  The  method  then  iteratively 
calculates  strong  state  gradients,  target  track  direction,  and 
updates  strong  state  vectors  over  the  analysis  grid. 

Probability  density  functions  are  recalculated.  Calculations 
are  reiterated  if  the  probability  density  functions  change  more 
than  an  information  threshold.  The  strong  state  vectors  and 
weak  state  probability  density  functions  are  provided  as  output 
if  the  information  threshold  is  not  exceeded. 

BRIEF  DESCRIPTION  OF  THE  DRAWINGS 
[0020]  Reference  is  made  to  the  accompanying  drawings  in 
which  are  shown  an  illustrative  embodiment  of  the  invention, 
wherein  corresponding  reference  characters  indicate 
corresponding  parts,  and  wherein: 

[0021]  FIG.  1  is  a  diagram  of  a  typical  target  motion 
analysis  system; 

[0022]  FIG.  2  is  a  graph  showing  target  motion  analysis 
parameters ; 

[0023]  FIG.  3  is  another  graph  showing  target  motion  analysis 
as  the  observer  and  target  move  through  time; 
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[0024]  FIG.  4.  is  a  graph  showing  the  cost  function  related 
to  the  initial  range  and  the  final  range; 

[0025]  FIG.  5  is  a  block  diagram  of  the  target  motion 
analysis  method; 

[0026]  FIG.  6A  and  6B  are  graphs  showing  target  motion 
analysis  under  the  prior  art  method;  and 

[0027]  FIG.  7A  and  7B  are  graphs  showing  target  motion 
analysis  under  the  current  method. 

DETAILED  DESCRIPTION  OF  THE  INVENTION 
[0028]  The  method  provided  herein  develops  a  more  robust 
target  motion  analysis  solution  estimate  and  provides  an 
improved  characterization  of  uncertainty  in  tracking  parameter 
estimates.  The  method  operates  over  a  tensor  product  space. 

Our  general  assumption  is  that  the  state  space  can  be 
partitioned  into  two  components  (or  subspaces) .  One  part  of  the 
state  space,  Ss,  corresponds  to  those  state  components  that  are 
strongly  observable.  Over  this  subspace,  parameter  estimates  can 
be  readily  obtained  via  gradient  based  estimation  methods  such 
as  steepest  descent,  Gauss-Newton  and  Newton-Raphson .  In  the 
other  part  of  the  state  space,  Sw,  the  state  vector  components 
are  considered  only  weakly  observable.  For  these,  applying 
gradient  based  estimation  methods  can  be  problematic  with  regard 
to  gradient  stability  and  algorithm  convergence.  The  method 
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herein  treats  the  strongly  observable  states  as  continuous 
random  variables  in  Ss  while  the  weakly  observable  states  are 
represented  by  a  collection  of  discrete  samples  distributed  over 
the  subspace  Sw.  The  tensor  product  space  comprises  the  cross 
product  of  these  continuous  and  sampled  variables  S  =  Ss  ®  Sw. 

The  method  herein  utilizes  a  non-linear  estimation  algorithm 
that  jointly  exploits  gradient  based  and  direct  search 
estimation  methods  in  a  manner  that  provides  a  robust 
algorithmic  stability  while  expediting  the  convergence  rate  of 
the  algorithm. 

[0029]  This  method  generates  an  efficient  and  informative 
uncertainty  representation  for  visualizing  solution  (state) 
sensitivity  over  the  weakly  observable  states.  A  depiction  of 
the  general  approach  of  the  invention  is  depicted  in  FIG.  5.  For 
general  application,  the  modeling  paradigm  might  include 
additional  parameters  whereby  unknown  modeling  parameters,  fn,  are 
adjoined  to  the  strongly  observable  state  components,  xs,  such 
that  for  each  weakly  observable  state  component  x^,,  k  =  1,  2...,  K, 
the  vector  xk  =  [xk,fn]  that  minimizes  the  performance  index  is 
determined.  Conditioning  on  the  weakly  observable  states 
improves  the  efficiency  of  the  estimation  process.  It  is  assumed 
that  any  modeling  parameters  fn  that  are  adjoined  to  the  strongly 
observable  states,  xs,  will  yield  an  estimation  vector  x  of 
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dimension  N  that  remains  strongly  observable  such  that  gradient 
based  estimation  methods  remain  applicable. 

[0030]  Parameter  estimates  of  xk  are  formed  for  each  weak 
state  sample  point  by  applying  the  non-linear  gradient  based 
estimation  methods  to  the  observation  data  Z.  This  provides  a 
stable  mechanism  for  maximizing  the  goodness-of-f it  to  the 
observed  data  at  each  of  the  sample  points.  As  a  consequence, 
parameter  estimates  of  xk  will  vary  across  the  set  of  sample 
points,  xk .  In  this  way,  the  plausibility  of  each  of  the  sample 
points  is  maximized.  These  estimates  are  then  used  to  develop  a 
measure  of  the  uncertainty  over  the  weakly  observable  subspace 
S„. 

[0031]  When  our  improved  estimation  method  is  applied  to  the 
endpoint  coordinate  system  without  any  additional  modeling 
parameters,  the  endpoint  bearings  constitute  the  strongly 
observable  states,  xs  =  W1,B2\r  as  tiedown  bearings  are  directly 
observable.  The  standard  approach  to  the  PEP  algorithm  is  to 
adopt  measured  bearings  at  the  prescribed  endpoint  times.  Target 
range  must  be  determined  from  the  bearing  data  and,  depending 
upon  the  geometry  applied  to  collect  the  data,  can  be  weakly 
observable.  Hence,  xw  =  [Ri,R2]  denotes  the  weakly  observable 
component  of  the  state  vector  that  is  sampled  over  the  Rlt  R2  grid 
associated  with  the  PEP.  When  the  bearing  tiedowns  are  optimized 


14  of  25 


Attorney  Docket  No.  300118 


for  each  sample  x^,,  a  distinct  set  of  bearing  tiedowns  may  be 
calculated  for  each  Rlt  R2  grid  point.  Our  approach  to  improved 
uncertainty  representation  in  this  endpoint  range  space  is  to 
apply  these  distinct  values  in  the  calculation  of  the  PEP  RMS 
cost  function  that  is  displayed. 

[0032]  The  approach  above  improves  on  the  typical  grid 
approach  to  displaying  uncertainty  in  that  refinement  in  the 
observed  parameters  gives  a  more  informative  slice  of  the  total 
state  space  uncertainty.  The  method  herein  gives  a  more  robust 
representation  of  the  range  parameters  than  in  the  case  of  the 
PEP  surface  shown  in  FIG.  4.  Once  the  subspace  likelihood  has 
been  expanded  to  its  plausible  extent  via  maximum  likelihood 
estimation  of  x,  the  parametric  uncertainty  can  be  projected  to 
any  other  subspace  of  interest. 

[0033]  Computational  efficiency  is  obtained  by  noting  that 
not  all  the  grid  points  yield  the  same  impact  on  the  uncertainty 
representation.  It  is  inefficient  to  iterate  extensively  in 
order  to  fit  the  best  observed  parameters  for  each  grid  point 
independent  of  their  contribution  to  the  uncertainty 
representation.  grid  points  that  correspond  to  implausible 

or  unlikely  tracks  will  have  little  impact  on  the  overall 
uncertainty  representation.  To  address  this,  each  step  in  the 
estimation  process  is  performed  simultaneously  for  each  sampled 
grid  point.  At  each  iteration  in  the  estimation  process,  the 
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likelihood  function  that  is  produced  over  the  weak  state 
subspace  Sw  given  the  current  iteration  is  calculated  and 
compared  to  that  attained  in  the  previous  iteration.  The 
quantitative  measure  employed  for  this  comparison  is  the 
Kullback-Leibler  divergence.  It  is  known  in  many  information 
theoretic  and  statistical  estimation  contexts.  If  the 
difference  in  the  measure  over  the  grid  is  significant,  the 

iterative  estimation  process  is  allowed  to  continue.  If  it  is 
not,  the  process  is  terminated  thereby  saving  many  iteration 
calculations  that  might  be  executed  at  sampled  values  with  a 
poor  goodness-of-f it  measure. 

[0034]  This  approach  has  two  by-products.  One  is  that  it 
provides  a  completely  parallel  process  of  iterations  over  the 
grid.  This  promotes  computational  savings  as  parallelization 
methods  become  more  mature  in  software  development  products. 
Secondly,  it  has  the  effect  of  drastically  reducing  the  total 
iterations  required  since  tracks  (grid  points)  which  do  not 
significantly  impact  the  uncertainty  representation  can  be 
stopped  prior  to  achieving  a  local  minimum.  That  is,  unnecessary 
computation  is  not  utilized  to  calculate  implausible  tracks. 
[0035]  This  method  incorporates  all  the  aspects  described 
above  with  some  very  important  extensions.  As  mentioned  above, 
the  state  vector  is  partitioned  into  two  components  governed  by 
how  strongly  or  weakly  the  component  vectors  are.  For  the  weakly 
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observable  components,  an  ensemble  of  samples  is  methodically 
drawn  over  the  weak  parameter  subspace.  With  augmentation  by  the 
strong  component  for  each  track,  an  ensemble  of  feasible  tracks 
is  developed.  The  strong  set  parameters  are  estimated  using  the 
standard  estimation  methods  to  maximize  the  plausibility  of  the 
weak  set.  In  developing  the  strong  set  parameters  using  gradient 
based  methods,  the  calculation  of  gradients  is  extended  from  a 
two  dimensional  array  over  the  measurement  set  and  the  states  of 
the  parameter  estimate  to  a  three  dimensional  matrix  formulation 
spanning  the  set  of  tracks  included  in  the  methodical  grid 
sampling  as  well  as  the  measurement  set  and  state  variable 
components . 

[0036]  In  this  extension.  A,  now  denotes  a  (K  *  M  *  N) 
gradient  matrix  of  partial  derivatives.  This  gradient  matrix  can 
be  formulated  such  that  each  row  corresponds  to  a  track  index  k, 
each  column  corresponds  to  a  measurement  index,  m.  Each  sheet 
is  associated  with  the  single  component,  x(n),  of  the  state 
variable  x.  Each  gradient  sheet  comprises  the  subset  of 
gradients , 


Xw) 


(17) 


3x(n) 

with  An  =  A(:,:,n).  The  track  parameter  index  n,  varies 
orthogonally  to  the  sheets. 
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Following  equation  15,  the  gradients  for  the  PEP  endpoint 
bearings  are  calculated  for  each  bearing  measurement  Bm  and  track 
combination  x£  as: 


where 


dB2 

VxO-Tri) 


(/?fcos(B1)Ar2ryfc(tm)  +  R,[sm(B1)AT2rk(trn))/rk(trn) 

(i?2  cos(B2)AT'1ry  (tm)  +  ^sin(B2)A7’1rx/c(tm))/r/c(tm) 
and  ry(tm)  are  calculated  as: 


(18) 


(19) 


r*  (tm)  =  A sin(B!)  +  A T2R\  sin(B2)  +  A x0(tm)  (20) 

ry{tm )  =  A T-lRI  cos(B!)  +  A T2R\  cos(B2)  +  Ay 0(tm)  (21) 

following  equations  11  and  12  with  range  at  time  tm  for  track  k 
given  by: 


rk(tm)  =  Jrk(tm)2  +  ryk(tm)2  ^ 

Note  that  the  observer  trajectory  deviations,  A x0(tm)  and  Ay0(tm) 
are  identical  for  each  track  k. 

[0037]  When  the  Householder  transformation  is  applied  to  this 
construct,  it  is  calculated  in  parallel  over  the  track  set 
{*w)/c=i  •  The  vectorised  Householder  is  applied  sequentially  over 
the  N  sheets  associated  with  x.  Following  equations  6  and  10  and 
using  a  step-size  of  unity  (s=l  in  equation  6) ,  the 
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transformation  operates  over  the  measurement  set  to  yield  the 
vectorised  back  substitution: 


gf  =  (Rk')~1zk  (23) 

for  k—l,...,K. 

[0038]  The  second  significant  extension  is  to  control  the 
iteration  process  over  which  the  method  converges  to  the  optimal 
solution  set.  To  do  this,  the  method  measures  the  change  in 
information  that  is  gained  over  the  weak  observability  parameter 
subspace,  Sw,  each  time  that  a  set  of  track  parameters  are 
updated.  During  each  iteration,  an  estimate,  xk ,  of  the 
parameters  in  xk  that  are  associated  with  track  k  are  updated 
according  to 


Xi+1=  xk  +  gk,k  =  (24) 

[0039]  The  resulting  likelihood  function  developed  over  the 
subspace  Sw  to  a  conditional  probability  density  function  f(k)  by 
assuming  a  diffuse  prior  over  Sw  and  normalizing  the  result  such 

that  Z»=i/t(fc)  =  i. 


/t(fc)  = 


,-M*n 


(25) 


The  Kullback-Leibler  divergence  is  then  calculated  between  the 
probability  density  function  produced  in  the  preceding 
iteration,  f,  and  its  update  following  the  iteration,  f+ 1. 
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[0040]  The  Kullback-Leibler  divergence  is  computed  as: 


D(fi(.k)\fi+1(k))  =  ^  fi(k)logK  (^^)) 

/c — 1 

When  the  change  in  information  falls  below  a  specified 
threshold. 


(26) 


<  e.  (27) 

the  iteration  is  halted.  The  information  measure  is  developed 
following  Shannon'’  s  practice  of  applying  the  logarithm  base 
equal  to  the  size  of  the  input  alphabet  to  avoid  scaling  issues. 
This  does  not  pose  difficulty  as  any  base  logarithm  can  be 
applied;  that  is  loga(p)  =  logb  (p) /logb  (a)  .  The  equivalent 
termination  criterion  using  a  logarithmic  base  of  2  becomes 

(a.  s  *  (28) 

D(fi(k)\fi+1(k))  =  log2  (-^-J  +  ao^/iCWO^+i)  “/(**))  <  elo92K 

1  k=  1 

where 


ai  =  £e-M*n 


fff+i 


k= 1 


(29) 


(30) 


a0  —  0.5 


log2(e ) 


er 


(31) 
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With  this  approach,  a  fixed  termination  criterion  can  be 
applied.  The  value  of  £=0.001. 

[0041]  A  functional  diagram  of  the  invention  is  provided  in 

FIG.  5.  This  diagram  illustrates  the  dichotomy  that  is  created 
between  the  strongly  and  weakly  observable  state  components.  In 
the  current  embodiment  bearings  are  strongly  observable  state 
components  and  ranges  are  weakly  observable  state  components; 
however,  other  strongly  and  weakly  observable  state  components 
can  be  used.  Signal  strength  could  be  an  additional  strongly 
observed  component.  An  example  of  an  additional  weakly 
observable  component  is  target  velocity. 

[0042]  In  FIG.  5  step  30  these  strongly  and  weakly  observable 
components  are  received  at  the  target  motion  analysis  system  24 
from  a  processor  such  as  20  of  FIG.  1.  In  step  32  the  weakly 
observable  components  Sw  are  utilized  to  form  an  analysis  grid 
that  can  be  used  for  parallelizing  computation.  Strongly 
observed  components  Ss  are  initialized  to  initial  conditions  in 
step  34.  In  step  36  the  initial  Ss  vectors,  Sw  analysis  grid, 
and  Sw  values  are  used  to  calculate  initial  Sw  probability 
density  functions. 

[0043]  Dashed  region  38  represents  computations  that  are 
performed  over  the  analysis  grid  formed  from  the  weakly 
observable  components  Sw  in  step  32 .  (Of  course  these 
computations  can  also  occur  sequentially.)  Calculation  of 
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strongly  observable  component  gradients  is  performed  in  step  40. 
This  step  utilizes  equation  (17)  discussed  above.  Step  42 
calculates  direction  of  the  strong  updates  utilizing  the 
vectorized  Householder  transformation  and  back-substitution  as 
provided  in  equation  (22) .  The  Ss  vectors  are  updated  in  step  44 
according  to  equation  (23)  . 

[0044]  After  completion  of  the  Sw  grid  calculations  of  38,  the 
method  determines  if  calculations  are  converging  to  an  optimal 
solution  set.  In  step  46,  the  Sw  probability  density  functions 
are  updated  in  accordance  with  equation  (24) .  Step  48 
calculates  the  Kullback-Liebler  divergence  between  the  Sw 
probability  density  functions  from  the  prior  step  and  those  of 
the  current  step  to  determine  the  change  in  information  over  the 
weakly  observable  data  subspace.  This  change  in  information  is 
compared  against  a  threshold  e  in  step  50.  (The  threshold  can 
be  either  calculated  or  pre-established  in  accordance  with  pre- 
established  tolerances  as  discussed  above  related  to  equation 
(27) -(31) .)  If  the  change  in  information  is  above  the  threshold 
e,  updated  parameters  are  provided  to  the  parallelized  routines 
38  in  step  52  for  additional  iterative  computation.  If  the 
change  in  information  falls  below  the  threshold  e  the  Ss  vectors 
and  Sw  likelihood  are  provided  as  output  in  step  54. 

[0045]  An  illustration  of  the  impact  of  the  invention  is 
provided  in  FIGS.  6A,  6B,  7A,  and  7B.  Here,  the  method  is 
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applied  to  perform  smoothing  of  the  endpoint  bearings  over  the 
estimated  Ri,  R2  grid  and  to  construct  the  area  of  uncertainty  in 
geographic  coordinates.  FIGS.  6A  and  7A  provide  a  geographic 
display  of  observer  and  target  tracks  for  a  bearings  only  target 
motion  analysis  example.  The  observer  trajectory  consists  of 
two  legs  identified  as  60A  and  60B  in  FIG.  6A  and  as  70A  and  70B 
in  FIG.  7A.  Endpoint  bearings  Bi  and  B2  emanate  from  the  terminal 
points  of  the  observer  trajectory.  The  true  target  track  is 
shown  as  62  in  FIG.  6A  and  as  72  in  FIG.  7A.  FIGS.  6B  and  7B 
provide  an  illustration  of  the  RMS  cost  surface  developed  as  a 
function  of  the  endpoint  ranges.  The  minimal  cost  region  66 
shows  up  as  the  long  slightly  tilted  region  comprised  by  the 
collection  of  'x's  at  the  center  of  the  fan  of  expanding  wedges. 
FIGS.  6A  and  6B  illustrate  the  result  obtained  prior  to 
application  of  the  method  where  bearings  are  tied  down 
arbitrarily.  The  minimum  RMS  solution  is  indicated  as  track  line 
64  that  broaches  the  endpoint  bearings  in  the  geometry  plot.  The 
corresponding  AOU  is  illustrated  as  the  set  of  squares  on  the 
second  bearing  line  B2 .  FIG.  7A  illustrates  the  impact  of  the 
method.  Once  the  current  method  has  calculated  the  endpoint 
bearings  at  each  track  point  (i.e.,  each  range  1  and  range  2 
node  within  the  range  grid) ,  the  set  of  points  in  each  of  the 
RMS  shading  bands  expands  to  provide  a  more  accurate  depiction 
of  the  uncertainty  in  the  range  estimates.  The  improvement  in 
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the  estimated  target  track  74  due  to  the  endpoint  smoothing  is 
evident  in  the  figure.  Also  note  that  each  successive 
sensitivity  level  has  increased  in  size.  The  inner  region 
identified  by  squares  comprising  the  AOU  is  shown  to  include  the 
true  target  track  whereas  in  FIG.  6A,  it  did  not.  In  FIG.  7B, 
the  minimal  cost  region  is  indicated  by  reference  number  76. 
[0046]  The  utilization  of  this  method  has  several  advantages 
over  traditional  data  processing  methods.  By  constructing  the 
track-measurement-parameter  multi-dimensional  array,  vectorized 
numerical  operations  can  be  applied  over  the  ensemble  of 
candidate  target  tracks  to  increase  the  efficiency  of  the 
algorithm.  By  decomposing  the  target  state  vector  into  strongly 
and  weakly  observable  components,  the  sensitivity  of  the 
objective  function  to  additional  modeling  parameters  can  be 
evaluated  by  adjoining  these  parameters  to  the  strong  state 
components  and  applying  gradient  based  estimation  methods.  By 
devising  an  estimation  termination  criteria  that  measures  the 
change  of  information  content  over  the  weakly  observable 
parameter  subspace,  parameter  estimation  is  allowed  to  converge 
only  for  those  track  candidates  that  actually  influence  the 
determination  of  the  uncertainty  region  boundary. 

[0047]  The  current  method  allows  a  more  robust  uncertainty 
representation  with  a  reasonable  computational  burden.  The 
novelty  here  is  in  the  fact  that  maximum  likelihood  estimation 
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is  applied  simultaneously  to  the  collection  of  feasible  tracks 
rather  than  individually  to  each  track  with  no  regard  to 
convergence  properties  for  the  track.  This  significantly 
expedites  the  convergence  of  the  overall  algorithm. 

[0048]  It  will  be  understood  that  many  additional  changes  in 
the  details,  materials,  steps  and  arrangement  of  parts,  which 
have  been  herein  described  and  illustrated  in  order  to  explain 
the  nature  of  the  invention,  may  be  made  by  those  skilled  in  the 
art  within  the  principle  and  scope  of  the  invention  as  expressed 
in  the  appended  claims. 

[0049]  The  foregoing  description  of  the  preferred  embodiments 
of  the  invention  has  been  presented  for  purposes  of  illustration 
and  description  only.  It  is  not  intended  to  be  exhaustive,  nor 
to  limit  the  invention  to  the  precise  form  disclosed;  and 
obviously,  many  modification  and  variations  are  possible  in 
light  of  the  above  teaching.  Such  modifications  and  variations 
that  may  be  apparent  to  a  person  skilled  in  the  art  are  intended 
to  be  included  within  the  scope  of  this  invention  as  defined  by 
the  accompanying  claims. 
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UNCERTAINTY  PREDICTION  IN 
PASSIVE  TARGET  MOTION  ANALYSIS 

ABSTRACT  OF  THE  DISCLOSURE 

A  method  for  target  motion  analysis  for  calculating  the 
track  between  an  observer  and  a  target  and  estimating  its 
uncertainty  includes  initializing  strongly  observable  state 
components  into  strong  state  vectors  and  weakly  observable  state 
components  to  form  a  weak  state  analysis  grid.  Initial 
probability  density  functions  are  calculated  from  all 
components.  The  method  then  iteratively  calculates  strong  state 
gradients,  target  track  direction,  and  updates  strong  state 
vectors  over  the  analysis  grid.  Probability  density  functions 
are  recalculated.  Calculations  are  reiterated  if  the 
probability  density  functions  change  more  than  an  information 
threshold.  The  strong  state  vectors  and  weak  state  probability 
density  functions  are  provided  as  output  if  the  information 


threshold  is  not  exceeded. 
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